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This paper deals with efficient numerical representation and manipulation of differential and 
integral operators as symbols in phase-space, i.e., functions of space x and frequency ^. The 
symbol smoothness conditions obeyed by many operators in connection to smooth linear partial 
differential equations allow to write fast-converging, non-asymptotic expansions in adequate sys- 
tems of rational Chebyshev functions or hierarchical splines. The classical results of closedness 
r^ of such symbol classes under multiplication, inversion and taking the square root translate into 

"t^ practical iterative algorithms for realizing these operations directly in the proposed expansions. 

Because symbol-based numerical methods handle operators and not functions, their complexity 
depends on the desired resolution N very weakly, typically only through log N factors. We 
present three applications to computational problems related to wave propagation: 1) precon- 
ditioning the Helmholtz equation, 2) decomposing wavefields into one-way components and 3) 
depth-stepping in reflection seismology. 
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A typical problem of interest in this paper is the efficient representation of functions of elliptic 
linear operators such as 
(^ A = /-div(a(a;)V-), 

where a{x) > c > is smooth, and x G [0,1]^ with periodic boundary conditions. We have in 
mind the inverse, the square root, and the exponential of A as important examples of functions of 
A. While most numerical methods for inverting A, say, would manipulate a right-hand side until 
convergence, and leverage sparsity or other properties of A in doing so, the scope of this paper is 
quite different. Indeed, we present expansions and iterative algorithms for manipulating operators 
as "symbols" , which make no or little reference to the functions of x to which these operators may 
later be applied. 

The central question is of course that of choosing a tractable representation for differential 
and integral operators. If a function f(x) has N degrees of freedom — if for instance it is sampled 
on N points — then a direct representation of operators acting on f{x) would in general require 
N'^ degrees of freedom. There are many known methods for bringing down this count to 0{N) or 
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0(A^log A^) in specific cases, such as leveraging sparsity, computing convolutions via FFT, low-rank 
approximations, fast summation methods |l19ji21j, wavelet or x-let expansions ^2j, partitioned SVD 
and H-matrices [U |20] . 

The framework presented in this paper is different, in the sense that we aim for a complexity 
essentially independent of A^, i.e., at most a low-degree polynomial of log A^, for representing and 
combining operators belonging to standard classes. Like the methods above, the operator is "com- 
pressed" in such a way that applying it to functions remains simple; it is only for this operation 
that the complexity needs to be greater than N, in our case O(A^logA^). 

1.1 Smooth symbols 

Let A denote a generic differential or singular integral operator, with kernel representation 

id 



Af{x) = / k{x,y)f{y)dy, x,y e 

Expanding the distributional kernel k{x, y) in some basis would be cumbersome because of the 
presence of a singularity along the diagonal x = y. For this reason we choose to consider operators 
as pseudodifferential symbols a(x,^), by considering their action on the Fourier transfornrJ/(^) of 
fix); 

Af{x) = j e^^'-Mx,i)f{i)di. 

One often writes a{x,D) for A, where D = —iVx/i^n). In this representation, the singularities of 
k{x,y) are turned into the oscillating factor g^'^*^''^, which can be discounted by focusing on the 
symbol a(x,^). The latter is usually smooth, and in a very peculiar way. It is the special form 
of the smoothness estimates for a — which we now describe — that guarantees the efficiency of the 
discretizations proposed in this paper. 

A symbol defined on M^ x W^ is said to be pseudodifferential of order m (and type (1,0)) if it 
obeys 

\d^d^a{x,o\ < c„;3(0™-i"i, (0 ^ (1 + ICP)'/', (1) 

for all multi-indices a,p. This symbol class is denoted 5""; the operator corresponding to some a 
in this class is denoted a(x, D), and belongs by definition to the class ^™'. Manifestly, one power of 
(0 is gained for each differentiation, meaning that the larger (^), the smoother a. For instance, the 
symbols of differential operators are polynomials in ^ and obey ([I]) when they have C°° coefficients. 
Large classes of singular integral operators also have symbols in the class 5"" |35] . 

The standard treatment of pseudodifferential operators makes the further assumption that some 
symbols can be represented as poly homogeneous series, such as 

a(x,e)~^a,(x,argOler-^ (2) 

j>o 

which defines the "classical" symbol class 5™ when the Uj are of class C°° . Corresponding operators 
are said to be in the class ^™. The series should be understood as an asymptotic expansion; it 
converges only when adequate cutoffs smoothly removing the origin multiply each term. Only then. 



^Our conventions in this paper: 

no = / e-'^'^-'fi^) dx. fix) = / e^^'^V(C) d^ 



the series does not converge to o(x, ^), but to an approximation that differs from a by a smoothing 
remainder r(x,^), smoothing in the sense that \dfdxr{x,C)\ = 0{{^)~°^). For instance, an operator 
is typically transposed, inverted, etc. modulo a smoothing remainder |24j . 

The subclass Q is central for applications — it is the cornerstone of theories such as geometrical 
optics — but the presence of remainders is a nonessential feature that should be avoided in the 
design of efficient numerical methods. The lack of convergence in Q may be acceptable in the 
course of a mathematical argument, but it takes great additional effort to turn such series into 
accurate numerical methods; see |3^ for an example. In a sense, the goal of this paper is to find 
adequate substitutes for ^ that promote asymptotic series into fast-converging expansions. 

There are in general no explicit formulas for the symbols of functions of an operator. Fortu- 
nately, some results in the literature guarantee exact closedness of the symbol classes ([I| or ([2]) 
under inversion and taking the square root, without smoothing remainders. A symbol a £ 5"", or 
an operator a{x, D) £ ^'", is said to be elliptic when there exists R > such that 

\a-\x,0\<C\^\-'^, when \C\ > R. 

• It is a basic result that if yl E *'"i, B € *™-2^ then AB G ^™-i+™2_ ggg fQj. instance Theorem 
18.1.8 in [21], Volumes. 



It is also a standard fact that if A G *™, then its adjoint A* G *'^. 

li A £ ^"*, and A is elliptic and invertiblqj on L^, then A~^ G ^"™. This result was proved 
by Shubinin 1978 in |33j. 



• For the square root, we also assume ellipticity and invertibility. It is furthermore convenient 
to consider operators on compact manifolds, in a natural way through Fourier transforms in 
each coordinate patch, so that they have discrete spectral expansions. A square root A^''^ 
of an elliptic operator A with spectral expansion A = ^- XjEj, where Ej are the spectral 
projectors, is simply 

A'/^ = Y,xfE„ (3) 

j 

with of course {A^''^)^ = A. In 1967, Seeley ^2] studied such expressions for elliptic A G *™, 
in the context of a much more general study of complex powers of elliptic operators. If in 
addition m is an even integer, and an adequate choice of branch cut is made in the complex 
plane, then Seeley showed that A^''^ G ^™ ; see [34J for an accessible proof that involves the 
complex contour "Dunford" integral reformulation of d3|. 

We do not know of a corresponding closedness result under taking the square root, for the non- 
classical class ^'". In practice, we will also manipulate operators that come from PDF on bounded 
domains with certain boundary conditions; the extension of the theory of pseudodifferential oper- 
ators to bounded domains is a difficult subject that this paper has no ambition of addressing. Let 
us also mention in passing that the exponential of an elliptic, non-self-adjoint pseudodifferential 
operator is not in general pseudodifferential itself. 

Numerically, it is easy to check that smoothness of symbols is remarkably robust under inversion 
and taking the square root of the corresponding operators, as the following simple one-dimensional 
example shows. 



^In the sense that A is a bijection from H"^{R'^) to L^(R''), hence obeys 11^/11^2 < C\\f\\H'^- EUipticity, in the 
sense in which it is defined for symbols, obviously does not imply invertibility. 



Let A := Att^I — div(a(x)V) on the periodic interval [0, 1] where a{x) is a random bandhniited 
function shown in Figure flla) . The symbol of this operator is 

a{x,0 = 47r2(l + a{x)\S,\'^) - 27rfVa(x) • ^, 

which is of order 2. In Figure II|b), we plot the values of a(x,^)(0~^ for x and ^ on a Cartesian 
grid. 

Since A is elliptic and invertible, its inverse C = A~'^ and square root D = A^'"^ are both 
well defined. Let use c(x, ^) and d{x, S.) to denote their symbols. From the above theorems, we 
know that the orders of c{x,^) and d{x,(,) are respectively —2 and 1. We do not believe explicit 
formulae exist for these symbols, but the numerical values of c(x, ^){C)'^ and d{x, £,){C)~^ are shown 
in Figure fife) and (d), respectively. These plots demonstrate regularity of these symbols in x and 
in ^; observe in particular the disproportionate smoothness in ^ for large |^|, as predicted by the 
class estimate ([I]). 
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Figure 1: Smoothness of the symbol in ^. (a) The coefficient a{x). (b) a(a;,^)(0~^ where a(x,^) 
is the symbol of A. (c) c(a;,^)(^)^ where c{x,^) is the symbol of C = A~^. (d) d{x,^){^)~^ where 
d{x, £,) is the symbol of -D = A^/'^. 



1.2 Symbol expansions 

Figure [T] suggests that symbols are not only smooth, but that they should be highly separable in x 
vs. ^. So we will use expansions of the form 

a(x,0 = 5^aA,^eA(x)5^(0(0'% (4) 

A 

where ex and g^ are to be determined, and (^)'^" = (1 + |^|2)<^a/2 encodes the order da of a{x,^). 
This choice is in line with recent observations of Beylkin and Mohlenkamp |1] that functions and 
kernels in high dimensions should be represented in separated form. In this paper we have chosen 
to focus on two-dimensional x, i.e. (x,^) G M^, which is already considered high-dimensional by 
the standards of numerical analysts. For all practical purposes the curse of dimensionality would 
prohibit any direct, even coarse sampling in M^. 

The functions ex{x) and 5^(^) should be chosen such that the interaction matrix ax^^ is as 
small as possible after accurate truncation. This choice also depends on the domain over which 
the operator is considered. In what follows we will assume that the rr-domain is the periodized 
unit square [0, 1]^ in two dimensions. Accordingly it makes sense to take for ex{x) the complex 
exponentials e^'^*^''*' of a Fourier series. The choice of gfj.{0 i^ more delicate, as x and ^ do not play 
symmetric roles in the estimate M. In a nutshell, we need adequate basis functions for smooth 
functions on M^ that behave like a polynomial of l/\(,\ as (, -^ cxd, and otherwise present smooth 
angular variations. We present two solutions in what follows: 

• A rational Chebyshev interpolant, where g^(^) are complex exponentials in angle 9 =arg ^, 
and scaled Chebyshev functions in |,^|, where the scaling is an algebraic map s = ^ ", . More 
details in Section [2. 1[ 



A hierarchical spline interpolant, where gfj.{C) are spline functions with control points placed 
in a multiscale way in the frequency plane, in such a way that they become geometrically 



scarcer as |^| — > oo. More details in Section 2.2 



Since we are considering x in the periodized square [0,1]^, the Fourier variable ^ is restricted 
to having integer values, i.e., ^ S Z^, and the Fourier transform should be replaced by a Fourier 
series. Pseudodifferential operators are then defined through 

a{x,D)fix) = Y, e2--«a(x,0/(6, (5) 

where /(^) are the Fourier series coefficients of /. That ^ is discrete in this formula should not be a 
distraction: it is the smoothness of the underlying function of ^ G M^ that dictates the convergence 
rate of the proposed expansions. 

The following results quantify the performance of the two approximants introduced above. We 
refer to an approximant as being truncated to M terms when all but at most M elements are put 
to zero in the interaction matrix a;^ ^ in d4| . 



Theorem 1. (Rational Chebyshev approximants). Assume that a £ 5™, that it is properly sup- 
ported, and assume furthermore that the aj in equation M) have tempered growth, in the sense 
that there exists Q, R > such that 

\d^d^aj{x,e)\<Q^^p-RK (6) 



Denote by a the rational Chebyshev expansion of a (introduced in Section \2.1^ , properly truncated 
to M terms. Call A and A the corresponding pseudodifferential operators on -^^([0, 1]^) defined by 
^. Then, there exists a choice of M obeying 

M <Cn- e"^/", Vn > 0, 
for some Cn > 0, such that 

\\A - ^||//m([o,l]2)^i2([o,l]2) < e. 

Theorem 2. (Hierarchical sphne approximants) . Assume that a € 5"*, and that it is properly 



supported. Denote by a the expansion of a in hierarchical splines for ^ (introduced in Section 2.2), 
and in a Fourier series for x, properly truncated to M terms. Call A and A the corresponding 
pseudodifferential operators on L^([0,1]^) defined by ([^. Introduce Pf<s the orthogonal projector 
onto frequencies obeying 

max(|Ci|,|6l)<A^. 

Then there exists a choice of M obeying 

M<C-e-2/{P+i).iog^r^ 

where p is the order of the spline interpolant, and for some C > 0, such that 

The important point of these theorems is that M is either constant in A^ (Theorem IT|, or grows 
Hke log N (Theorem^, where N is the bandlimit of the functions to which the operator is applied. 

1.3 Symbol operations 

At the level of kernels, composition of operators is a simple matrix-matrix multiplication. This 
property is lost when considering symbols, but composition remains simple enough that the gains 
in dealing with small interaction matrices a\^^ as in (Hi) are far from being offset. 

The twisted product of two symbols a and 6, is the symbol of their composition. It is defined as 
(att6)(x, D) = a(x, D)b{x, D) and obeys 

aS6(x, = j j e-^^'^^-yy^^-^^a{x, 7j)b{y, dydr^. 

This formula holds for £,,ri ^ M , but in the case when frequency space is discrete, the integral in 
r/ is to be replaced by a sum. In Section [3] we explain how to evaluate this formula very efficiently 
using the symbol expansions discussed earlier. 

Textbooks on pseudodifferential calculus also describe asymptotic expansions of aj^ft where neg- 
ative powers of 1^1 are matched at infinity |24| [T8| IM] . but, as alluded to previously, we are not 
interested in making simplifications of this kind. 

Composition can be regarded as a building block for performing many other operations using 
iterative methods. Functions of operators can be computed by substituting the twisted product for 
the matrix-matrix product in any algorithm that computes the corresponding function of a matrix. 
For instance, 

• The inverse of a positive-definite operator can be obtained via a Neumann iteration, or via a 
Schulz iteration; 



• There exist many choices of iterations for computing the square root and the inverse square 
root of a matrix [23^ , such as the Schulz-Higham iteration; 

• The exponential of a matrix can be obtained by the scahng-and-squaring method; etc. 

These examples are discussed in detail in Section [3| 

Two other operations that resemble composition from the algorithmic viewpoint, are 1) trans- 
position, and 2) the Moyal transform for passing to the Weyl symbol. They are also discussed 
below. 

Lastly, this work would be incomplete without a routine for applying a pseudodifferential oper- 
ator to a function, from the knowledge of its symbol. The type of separated expansion considered 
in equation (El) suggests a very simple algorithm for this task, detailed in Section ^ (This part is 
not original; it was already considered in previous work by Emmanuel Candes and the authors in 
|1U| . where the more general case of Fourier integral operators was considered.) 

1.4 Applications 

Applications of discrete symbol calculus abound in the numerical solutions of linear partial differ- 
ential equations (PDE) with variable coefficients. We outline several examples in this section and 
their numerical results are given in Section [4| 

In all of these applications, our solution takes two steps. First, we use discrete symbol calculus 
to construct the symbol of the operator which solves the PDE problem. Since the data has not 
been queried yet (i.e., the right hand side, the initial conditions, or the boundary conditions), the 
computational cost of this step is mostly independent of the size of the data. Once the operator is 
ready in its symbol form, we apply the operator to the data in the second step. 

The two regimes in which this approach could be preferred is when either 1) the complexity of 
the medium is low compared to the complexity of the data, or 2) the problem needs to be solved 
several times and benefits from being "preconditioned" in some way. 

A first, toy application of discrete symbol calculus is to the numerical solution of the simplest 
eUiptic PDE, 

Au := {I - div(a(a;)V)n = / (7) 

with a{x) > 0, and periodic boundary conditions on a square. If a{x) is a constant function, the 
solution requires only two Fourier transforms, since the operator is diagonalized by the Fourier 
basis. For variable a{x), discrete symbol calculus can be seen as a natural generalization of this 
fragile Fourier diagonalization property: we construct the symbol of A^^ directly, and once the 
symbol of A~^ is ready, applying it to the function / requires only a small number of Fourier 
transforms. 

The second application of discrete symbol calculus is related to the Helmholtz equation 

where the sound speed c(x) is a smooth function in x, in a periodized square. The numerical solution 
of this problem is difficult since the operator L is not positive definite so that efficient techniques 
such as multigrid cannot be used directly for this problem. A standard iterative algorithm, such 
as MINRES or BIGGSTAB, can easily take tens of thousands of iterations to converge. One way 
to obtain faster convergence is to solve a preconditioned system 

M-'^Lu = M-^f (9) 



with 



a;2 ^2 



M:=-A+^-- or M := -A + {1 + i)- 



Now at each iteration of the preconditioned system, we need to invert a hnear system for the 
preconditioner M. Multigrid is typicahy used for this [T3|, but discrete symbol calculus offers a 
way to directly precompute the symbol of M^^ . Once it is ready, applying M~^ to a function at 
each iteration is reduced to a small number of Fourier transforms — three or four when c{x) is very 
smooth — which we anticipate to be very competitive vs. a multigrid method. 

Another important application of the discrete symbol calculus is to "polarizing" the initial 
condition of a linear hyperbolic system. Let us consider the following variable coefficient wave 
equation on the periodic domain x € [0, 1]^, 

utt — div(a(a;)Vu) = 

u(0,x)=tio(x) (10) 

ut{^,x) = ui{x) 

with the extra condition J ui{x)dx = 0. The operator L := — div(a(2;)V) is symmetric positive 
definite, and let us define P to be its square root L^'"^. We can then use P to factorize the wave 
equation as 

{dt + iP){dt-iP)u = 0. 

As a result, the solution u{t,x) can be represented as 

u{t,x) = e^*^n+(x) + e-^*^n„(x) 

where the polarized components ti+(x) and u-{x) of the initial condition are given by 

uo + {iP)~^ui uo-{iP)~^ui 

ii_i = and U- = . 

+ 2 2 

To compute u+ and n_, we first use discrete symbol calculus to construct the symbol of P^^. Once 
the symbol of P~^ is ready, the computation of u+ and u_ requires only applying P^^ to the initial 
condition. Applying e**^ is a more difficult problem that we do not address in this paper. 

Finally, discrete symbol calculus has a natural application to the problem of depth extrapolation, 
or migration, of seismic data. In the Helmholtz equation 

A± + ^-^ + ^7 rU = F(x, z, k), 

OZ'^ C''[X,Z) 

we can separate the Laplacian as A = A_|_ + #^, and factor the equation as 

^ - Biz)^ V = F(x, z, k) - ^iz)u, (^ + B{z)^ u = v (11) 

where B = y^— A_|_ — ui"^ /c^{x, z) is called the one-way wave propagator, or single square root (SSR) 
propagator. We may then focus on the equation for f , called the SSR equation, and solve it for 
decreasing z from z = 0. The term ^{z)u above is sometimes neglected, as we do in the sequel, 
on the basis that it introduces no new singularities. 

The symbol of B^ is not elliptic; its zero level set presents a well-known issue with this type 
of formulation. In Section |4j we introduce an adequate "directional" cutoff strategy to remove the 



singularities that would otherwise appear, hence neglect turning rays and evanescent waves, and 
then use discrete symbol calculus to compute a well-behaved operator square root. We then show 
how to solve the SSR equation approximately using an operator exponential of B, also realized via 



discrete symbol calculus. Unlike traditional methods of seismic imaging (discussed in Section 1.6 
below), the only simplification we make here is the directional cutoff just mentioned. 



1.5 Harmonic analysis of symbols 

It is instructive to compare the symbol expansions of this paper with another type of expansion 
thought to be efficient for smooth differential and integral operators, namely wavelets. 

Consider x G [0, 1] for simplicity. The standard matrix of an operator ^ in a basis of wavelets 
'4'j,k{x) = 2^''^ijj{2^x — n) of L-^([0, 1]) is simply {ipj^k, Aipj/^f^'). Such wavelet matrices were first 
considered by Meyer in |29], and later by Beylkin, Coifman, and Rokhlin in [2,, for the purpose of 
obtaining sparse expansions of singular integral operators in the Calderon-Zygmund class. Their 
result is that either 0{N) or 0(A^log A^) elements suffice to represent a N-hy-N matrix accurately, 
in the ^2 sense, in a wavelet basis. This result is not necessarily true in other bases such as Fourier 
series or local cosines, and became the basis for much activity in some numerical analysis circles in 
the 1990s. 

In contrast, the expansions proposed in this paper assume a class of operators with symbols in 
the S"^ class defined (II|, but achieve accurate compression with 0(1) or O(logA^) elements, way 
sublinear in N. This stark difference is illustrated in Figure [2| 
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Figure 2: Left: the standard 512-by-512 wavelet matrix of the differential operator considered in 
Figured] truncated to elements greater than 10~^ (white). Right: the 65-by-15 interaction matrix 
of DSC, for the same operator and a comparable accuracy, using a hierarchical splines expansion in 
^. The scale is the same for both pictures. Notice that the DSC matrix can be further compressed 
by a singular value decomposition, and in this example has numerical rank equal to 3, for a singular 
value cutoff at 10~^. For values of N greater than 512, the wavelet matrix would increase in size 
in a manner directly proportional to A^, while the DSC matrix would grow in size like log N. 

Tasks such as inversion and computing the square root are realized in 0(log A^) operations, 
still way sublinear in A'^. It is only when the operator needs to be applied to functions defined on 
A^ points, as a "post-computation", that the complexity becomes C ■ NlogN. This constant C is 
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proportional to the numerical rank of the symbol, and reflects the difficulty of storing it accurately, 
not the difficulty of computing it. In practice, we have found that typical values of C are still much 
smaller than the constants that arise in wavelet analysis, which are often plagued by a curse of 
dimensionality |12j . 

Wavelet matrices can sometimes be reduced in size to a mere 0(1) too, with controlled accuracy. 
To our knowledge this observation has not been reported in the literature yet, and goes to show 
that some care ought to be exercised before calling a method "optimal" . The particular smoothness 
properties of symbols that we leverage for their expansion is also hidden in the wavelet matrix, as 
additional smoothness along the shifted diagonals. The following result is elementary and we give 
it without proof. 

Theorem 3. Let A G ^^ as defined by (Mj, for x G M and ^ G M. Let tpj^k be an orthonormal 
wavelet basis of L'^(R) of class C°° , and with an infinite number of vanishing moments. Then for 
each j, and each Ak = k — k' , there exists a function fj^Ak £ C°°(M) with smoothness constants 
independent of j, such that 

{'>Pj,k, Ail)j,k') = fjAki.'^'^k). 

We would like to mention that similar ideas of smoothness along the diagonal have appeared in 
the context of seismic imaging, for the diagonal fitting of the so-called normal operator in a curvelet 
frame |22| |9] . In addition, the construction of second-generation bandlets for image processing is 
based on a similar phenomenon of smoothness along edges for the unitary recombination of MRA 
wavelet coefficients [31] . We believe that this last "alpertization" step could be of great interest in 
numerical analysis. 

Theorem [3] hinges on the assumption of symbols in S™, which is not met in the more gen- 
eral context of Calderon-Zygmund operators (CZO) considered by Meyer, Beylkin, Coifman, and 
Rokhlin. The class of CZO has been likened to a limited-smoothness equivalent to symbols of type 
(1, 1) and order 0, i.e., symbols that obey 



\d?dSa{x,0\<CaAO' 



M+m 



Symbols of type (1,0) and order obeying ([I]) are a special case of this. Wavelet matrices of 
operators in the (1,1) class are almost diagonaPl but there is no smoothness along the shifted 
diagonals as in Theorem pj So while the result in [2] is sharp, namely no much else than wavelet 
sparsity can be expected for CZO, we may question whether the generality of the CZO class is 
truly needed for applications to partial differential equations. The authors are unaware of a PDE 
setup which requires the introduction of symbols in the (1,1) class that would not also belong to 
the (1,0) class. 

1.6 Related work 

The idea of writing pseudodifferential symbols in separated form to formulate various one-way 
approximations to the variable-coefficients Helmholtz equation has long been a tradition in seismic 
imaging. This almost invariably involves a high-frequency approximation of some kind. Some 



^Their standard wavelet matrix has at most 0{j) large elements per row and column at scale j — or frequency 
0(2-') — after which the matrix elements decay sufficiently fast below a preset threshold. L^ boundedness would follow 
if there were 0(1) large elements per row and column, but 0{j) does not suffice for that, which testifies to the fact 
that operators of type (1, 1) are not in general L^ bounded. The reason for this 0{j) number is that an operator 
with a (1, 1) symbol does not preserve vanishing moments of a wavelet — not even approximately. Such operators may 
turn an oscillatory wavelet at any scale j into a non-oscillating bump, which then requires wavelets at all the coarser 
scales for its expansion. 
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influential work includes the phase screen method by Fisk and McCartor pTj, and the generalized 
screen expansion of Le Rousseau and de Hoop t27j . This last reference discusses fast application 
of pseudodifFerential operators in separated form using the FFT, and it is likely not the only 
reference to make this simple observation. A modern treatment of leading-order pseudodifFerential 
approximations to one-way wave equations is in |37] . 

Expansions of principal symbols ao{x,^/\(,\) (homogeneous of degree is ^) in spherical har- 
monics in ^ is a useful tool in the theory of pseudodifFerential operators [SH] , and has also been used 
for fast computations by Bao and Symes in [1]. For computation of pseudodifFerential operators, 
see also the work by Lamoureux, Margrave, and Gibson |I26|. 

In the numerical analysis community, separation of operator kernels and other high-dimensional 
functions is becoming an important topic. Beylkin and Mohlenkamp proposed an alternated least- 
squares algorithm for computing separated expansions of tensors in [3l H], propose to compute 
functions of operators in this representation, and apply these ideas to solving the multiparticle 
Schrodinger equation in [5], with Perez. 

A different, competing approach to compressing operators is the "partitioned separated" method 
that consists in isolating ofF-diagonal squares of the kernel K(x, y) , and approximating each of them 
by a low-rank matrix. This also calls for an adapted notion of calculus, e.g., for composing and 
inverting operators. The first reference to this algorithmic framework is probably the partitioned 
SVD method described in [25]. More recently, these ideas have been extensively developed under 



the name H-matrix, for hierarchical matrix; see [SI ED] and http://www.hlib.org 

Separation ideas, with an adapted notion of operator calculus, have also been suggested for 
solving the wave equation; two examples are [6j and [13j. 

Exact operator square-roots — up to numerical errors — have in some contexts already been con- 
sidered in the literature. See [16^ for an example of Helmholtz operator with a quadratic profile, 
and [28] for a spectral approach that leverages sparsity, also for the Helmholtz operator. 

2 Discrete Symbol Calculus: Representations 

The two central questions of discrete symbol calculus are: 

• Given an operator A, how to represent its symbol a{x, S,) efficiently? 

• How to perform the basic operations of the pseudodifFerential symbol calculus based on this 
representation? These operations include sum, product, adjoint, inversion, square root, in- 
verse square root, and, in some cases, the exponential. 

These two questions are mostly disjoint; we answer the first question in this section, and the 
second question in Section |3] 

Let us write expansions of the form M|. Since e\{x) = e^'^*^ with x G [0, 1]^, we denote the 
x-Fourier series coefficients of a{x, ^) as 

axiO = I e-2^^^-^a(x,e) dx, A G I?. 

i[0,l]2 

We find it convenient to write ha^\{^) = oa(0(0~'^") hence 

a{x,0 = Y.^x[x)K,x{0{if^- (12) 
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In the case when a{-,^) is bandhniited with band Bx, i.e., ax{^) is supported inside the square 
{—Bx,Bx)'^ in the A-frequency domain, then the integral can be computed exactly by a uniform 
quadrature on the points Xp = p/{2Bx), with < pi,p2 < 2-6^. This grid is called X in the sequel. 

The problem is now reduced to finding adequate expansions ha^x for ha^X: either valid in the 
whole plane ^ G M^, or in a large square ^ E [— A^, A^]^. 

2.1 Rational Chebyshev interpolant 



For symbols in the class O), each function ha^xiO — ^a(0(0 '^^ is smooth in angle arg ^, and 
polyhomogeneous in radius |^|. This means that ha^x is for |^| large a polynomial of 1/|^| along 
each radial line through the origin, and is otherwise smooth (except possibly near the origin). 

One idea for efficiently expanding such functions is to map the half line |^| G [0, oo) to the 
interval [—1,1] by a rational function, and expand the result in Chebyshev polynomials. Put 
^ = {6,r), and fi = {m,n). Let 

where TLn are the rational Chebyshev functions [7 , defined from the Chebyshev polynomials of 
the first kind T„ as 

TLn{r) = T^{A-^\r)), 

by means of the algebraic map 

s I—, r = Al[s) = L- , r I— > s = Aj^ [r) = -. 

1 — S T -\- Ij 

The parameter L is typically on the order of 1. The proposed expansion then takes the form 

ha,x{0 = ^ax,f^g,,{0, 

or ha^xiO if the sum is truncated, where 

1 /■' /■'" -imer^ , ^ deds 



«A,M = ;^ / / Va((^, AL{s)))e-^"''Tn{s) 



For properly bandlimited functions, such integrals can be evaluated exactly using the right 
quadrature points: uniform in ^ £ [0,27r], and Gauss points in s. The corresponding points in 
r are the image of the Gauss points under the algebraic map. The resulting grid in the ^ plane 
can be described as follows. Let q = {qe,(lr) be a couple of integers such that ^ < qe < Ng and 
< qr < Nj.; we have in polar coordinates 

,„ .. ;2{AL{qr)-l] 

?g = 27r-— ,-cos 



We call this grid {^q} = il. Passing from the values ha^xiCq) to a^,^ and vice- versa can be done 
using the fast Fourier transform. Of course, ha^xiO is nothing but an interpolant of /ia,A(0 ^^ the 
points S,q- 

In the remainder of this section, we present the proof of Theorem [T] which contains the con- 
vergence rates of the truncated sums over A and /x. The argument hinges on the following L^ 
boundedness result, which is a simple modification of standard results in M'', see j35]. It is not 
necessary to restrict d = 2 for this lemma. 
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Lemma 1. Let a{x,S,) G C^ ([0, l]"',^oo(^'^)); where d' = d + 1 if d is odd, or d + 2 if d is even. 
Then the operator A defined by ([^ extends to a bounded operator on L'^{[0, l]'^), with 

The proof of this lemma is in the Appendix. 

Proof of Theorem^ Consider s = A~^ (r) G [—1,1) where A^ and its inverse were defined above. 
Expanding a{x,{6,r)) in rational Chebyshev functions TLn{r) is equivalent to expanding f{s) = 
a{x,{6.,Ai[s))) in Chebyshev polynomials Tn[s). Obviously, 

/oyl-ieC-([0,oo)) ^ /gC-([-l,l)). 

It is furthermore assumed that a(x, ^) is in the classical class with tempered growth of the 
polyhomogeneous components; this condition implies that the smoothness constants of /(s) = 
a{x, {6, Al{s))) are uniform as s — > 1, i.e., for all n > 0, 

3C„ : |/(")(5)I<C7„, se[-l,l], 

or simply, / G C°°([— 1, 1]). In order to see why that is the case, consider a cutoff function x(?') 
equal to 1 for r > 2, zero for < r < 1, and C°° increasing in between. Traditionally, the meaning 
of ([2]) is that there exists a sequence Sj > 0, defining cutoffs xi^^^j) such that 

a(x, (r, 6)) - Y, ajix, e)r-^x{rej) G S'^, 
i>o 

for all ?n, > 0. A remainder in 5"^°° = Um>o '^cl"* ^^ called smoothing. As long as the choice of 
cutoffs ensures convergence, the determination of a(x,^) modulo S~°^ does not depend on this 
choice. (Indeed, if there existed an order-fc discrepancy between the sums with x(rej) or x(^^i)) 
with k finite, it would need to come from some of the terms ajr~^{x{r£j) — x(r(5j)) for j < k. But 
each of these terms is of order —00.) 

Because of condition (|6J), it is easy to check that the particular choice Ej = 1/{2R) suffices for 
convergence of the sum over j to a symbol in S^. As mentioned above, changing the Ej only affects 
the smoothing remainder, so we may focus on Ej = 1/{2R). 

After changing variables, we get 

/(.) = a{x, {e, AUs))) = j; a,(x, 6)1-^ H^Yx (^) + r{s), 

j>o V "T / V / 

where the smoothing remainder r{s) obeys 

|r(")(s)|<C7„,M(l-s)^^ VM>0, 

hence, in particular when M = 0, has uniform smoothness constants as s ^ 1. It suffices therefore 
to show that the sum over j > can be rewritten as a Taylor expansion for /(s) — r(s), convergent 
in some neighborhood of s = 1. 

Let z = 1 — s. Without loss of generality, assume that R > 2L, otherwise increase R to 2L. 
The cutoff factor x ( ^2R ) equals 1 as long as < 2: < j^. In that range. 



f{\-z)-r{\-z) = Y,a,{x,Q)L"'^ 



yl 



(2 - z)i 
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By making use of the binomial expansion 

(2 - z)i ~ 2^ UJ V i - 1 

^ ' m>0 ^ 

and the new index k = j + m, we obtain the Taylor expansion about z = 0: 






k>0 l<j<k 



^-1\ < v-fc-i (k- i\ _ ^k-1 
j -ij - ^"=0 I n 



To check convergence, notice that ( . -, ) < Sn=o I ^ ) = 2 , combine this with (6), and 
obtain 



2-k Y^ aj{x,e) //c - 1\ Qoo y^ f R^ ^ 

l<j<k ^ ^ l<j<A: ^ 



2 1- L/R \L 

We assumed earlier that z G [0,L/(4i?)]: this condition manifestly suffices for convergence of the 
sum over k. This shows that / E C°°{[—1, 1]); the very same reasoning with Q^/s in place of Qoo 
also shows that any derivative d^dg f{s) G C°°([— 1, 1]). 

The Chebyshev expansion of /(s) is the Fourier-cosine series of /(cos0), with (j) £ [0,7r]; the 
previous reasoning shows that /(cos0) G C°°([0, cxd]). The same is true for any {x,9) derivatives 
of /(cos(/)). 

Hence a(x, (ylL(cos (j)),6)) is a C°° function, periodic in all its variables. The proposed expansion 
scheme is simply: 

• A Fourier series in x G [0, 1]^; 

• A Fourier series in ^ G [0,27r]; 

• A Fourier-cosine series in (p G [0,7r]. 

An approximant with at most M terms can then be defined by keeping [M^'^J Fourier coefficients 
per direction. It is well-known that Fourier and Fourier-cosine series of a C°°, periodic function 
converge super-algebraically in the L°° norm, and that the same is true for any derivative of the 
function as well. Therefore if om is this M-term approximant, we have 

sup|af(a-a)(x,(ylL(cos0),0))| <Cp,M-M-'^, V multi-index /3 

x,6,<j> 

We now invoke Lemma 111 with a — um in place of a, choose M = 0{e^^''^) with the right 
constants, and conclude. 

D 

It is interesting to observe what goes wrong when condition (Im is not satisfied. For instance, if 
the growth of the aj is fast enough in Q, then it may be possible to choose the cutoffs x(^jlCl) such 
that the sum over j replicates a fractional negative power of \^\, like |^|~^'^, and in such a way that 
the resulting symbol is still in the class defined by (IT]) . A symbol with this kind of decay at infinity 
would not be mapped onto a C°° function of s inside [—1,1] by the algebraic change of variables 
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Al, and the Chebyshev expansion in s would not converge spectrally. This kind of pathology is 
generally avoided in the literature on pseudodifferential operators by assuming that the order of 
the compound symbol a{x,S,) is the same as that of the principal symbol, i.e., the leading-order 
contribution ao(a;,arg ^). 

Finally, note that the obvious generalization of the complex exponentials in arg ^ to higher- 
dimensional settings would be spherical harmonics, as advocated in ^IJ. The radial expansion 
scheme should probably remain unchanged, though. 

2.2 Hierarchical spline interpolant 

An alternative representation is to use a hierarchical spline construction in the ^ plane. We define 
ha,x{0 to be an interpolant of ha,x{C) = o,\{0{0~'^° ^ follows. We only define the interpolant in 
the square (, G [— A^, A^]^ for some large A. Pick a number B^ — independent of A — that plays the 
role of coarse-scale bandwidth; in practice it is taken comparable to B^. 



• Define Do = {-B^,B^)\ For each ^ G Dq, KAO ■= KAO- 

• For each j = 1,2,--- ,L = log^^iN/B^), define Dj = {-3^B^,3^B^f - Dj^i. 



We further 



partition Dj into eight blocks: 



D, 



U ^^V 



i=l 



where each block Dj^i is of size 2 • 3^^^B^ x 2 • 3^^^B^. Within each block Dj^i, we sample 
ha,x{0 with a Cartesian grid Gj^i of a fixed size. The restriction of ha,x{0 i^ ^j,i is defined 
to be the spline interpolant of ha^x{S,) on the grid Gj,j. 
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(b) 



Figure 3: Hierarchical spline construction. Here B^ = Q, L = 4, and A^ = 486. The grid Gj^i is of 
size 4x4. The grid points are shown with "-I-" sign, (a) The whole grid, (b) The center of the 
grid. 

We emphasize that the number of samples used in each grid Gj^i is fixed independent of the 
level j. The reason for this is that the function ha^xiO g^-ins smoothness as S, grows to infinity. In 
practice, we set Gj^i to be a 4 x 4 or 5 x 5 Cartesian grid and use cubic spline interpolation. 
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Let us summarize the construction of the representation a(x,^) ss J2\^>^i^)^a,x{0{0'^°'- ^s 
before fix a parameter Bx that governs the bandwidth in x, and define 






,P2 < 2Bx\ and ^ = Do\J l[j Gj,i j 



The construction of the expansion of a(x, ^) takes the following steps 

• Sample a{x, ^) for all pairs of {x, ^) with x £ X and ^ G $7. 

• For a fixed ^ G $7, use the fast Fourier transform to compute ax{S.) for all A G {—Bx, Bx)'^- 



• 



For each A, construct the interpolant /ia,A(0 from the values of ha,x{£,) = aA(0(0 ^-t 
e G S7. 



Let us study the complexity of this construction procedure. The number of samples in X 
is bounded by 4i?^, considered a constant with respect to N. As we use a constant number of 
samples for each level j = 1,2,- ■ ■ ,L = log^{N/B^), the number of samples in J7 is of order 
O(logA^). Therefore, the total number of samples is still of order O(logA^). Similarly, since the 
construction of a fixed size spline interpolant requires only a fixed number of steps, the construction 
of the interpolants {/ia,A(0} takes only 0(log A^) steps as well. Finally, we would like to remark 
that, due to the locality of the spline, the evaluation of /ia,A(0 for any fixed A and (, requires only 
a constant number of steps. 

We now expand on the convergence properties of the spline interpolant. 

Proof of Theorem [^ If the number of control points per square Dj^i is K^ instead of 16 or 25 as 
we advocated above, the spline interpolant becomes arbitrarily accurate. The spacing between two 
control points at level j is 0(3^ /K). With p be the order of the spline scheme — we took p = 3 
earlier — it is standard polynomial interpolation theory that 

sup \ha,\{0 - K,x{0\ < Ca,\.p • T7 • SUp ||(9?/la,A||L°o(D,- ,)• 



The symbol estimate M guarantees that the last factor is bounded by C ■ suptg^). .(0 ^ ^. Each 
square Dj^i, for fixed j, is at a distance 0(3^) from the origin, hence sup^^^ . .{^)~P'^ = 0(3^-'^^+-'^)). 
This results in 

sup \ha,x{0 - ha,xm < Ca,X,p ' K'^'^ . 



^GD 



3,1 



This estimate is uniform over -Dj,i, hence also over ^ G [— A, A]^. As argued earlier, it is achieved 
by using O(Ar^logA^) spline control points. If we factor in the error of expanding the symbol in 
the X variable using AB"^ spatial points, for a total of M = 0{B'^K'^ log N) points, we have the 
compound estimate 

sup sup \a{x,0-a{x,0\<C-{B~'^ + K-P~^). 
The same estimate holds for the partial derivatives of a — a in x. 
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Functions to which the operator defined by a{x, ^) is apphed need to be bandhmited to [—N, iV]^, 
i.e., /(.^) = for ^ ^ [— A^, A^]^, or better yet / = P/v/. For those functions, the symbol a can be 
extended by a outside of [— A^, -/V]^, Lemma 111 can be apphed to the difference A — A, and we obtain 

\\iA-A)fU2<C-{B-^+K-P-').\\fh2 

The leading factors of ||/||l2 in the right-hand side can be made less than e if we choose B = 
Q^^-i/oo^ and K = 0{e~^'^^'^^'), with adequate constants. The corresponding number of points 
in X and ^ is therefore M = 0(e^^/(P+^) • log TV). 

D 

3 Discrete Symbol Calculus: Operations 

Let A and B be two operators with symbols a{x,^) and b{x,S,)- Suppose that we have already 
generated their expansions 

a{x,0 « j;eA(x)/ia,A(e)(0''" and b{x,0 « J] eA(x)/lb,,(0(0''^ 

A A 

where da and db are the orders of a{x,^) and b{x,S,), respectively. It is understood that the sum 
over A is truncated, that ha,\{^) are approximated with /iq.aC?) by either method described earlier, 
and that we will not keep track of which particular method is used in the notations. 
Let us now consider the basic operations of the calculus of discrete symbols. 

Scaling C = aA. For the symbols, we have c{x, ^) = aa{x, ^). In terms of the Fourier coefficients, 

CA(e)=aaA(e)^«VA(0(0'"- 

Then dc = da and 

KAO ■■= haAO- 

Sum C = A + B. For the symbols, we have c(x,.^) = a(x,^) + b{x,^). In terms of the Fourier 
coefficients, 

£a(0 = ax{0 + bxiO « Va(0(0'" + hAOiO"'- 

Then dc = niax{da, d^) and h^AC) is the interpolant that takes the values 



at ^ G i7, where $7 is either the Gaussian points grid, or the hierarchical spline grid defined earlier. 
Product C = AB. For the symbols, we have 

r, J 

In terms of the Fourier coefficients, 

k+l=X k+l=X 
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Then dc = da + db and /ic,A(0 is the interpolant that takes the values 

\fc+/=A / 

at ^ G O. 

Transpose C = A* . For the symbols, it is straightforward to derive the formula 
In terms of the Fourier coefficients, 



cxiO = a_A(e + A) « KMi + A)(e + A)'^''. 
Then dc = da and hc\{S,) is the interpolant that takes the values 

at ^ G il. 

Inverse C = j4~^ where ^4 is symmetric positive definite. We first pick a constant a such that 
a\a{x,^)\ ^ 1 for ^ G {—N,N)'^. Since the order of a(x,^) is da, a ss ©(l/A^*^"). In the following 
iteration, we first invert aA and then scale the result by a to get C. 

• Xo = I. 

• For k = 0,1,2, .. ., repeat X/^^i = 2Xfc — ^^(q^)^^ until convergence. 

• Set C = aXk- 

This iteration is called the Schulz iteration, and is quoted in [3 . It can be seen as a modified 
Newton iteration for finding the nontrivial zero of f{X) = XAX — X, where the gradient of / is 
approximated by the identity. 

As this algorithm only utilizes the addition and the product of the operators, all of the compu- 
tation can be carried out via discrete symbol calculus. Since a ~ 0{1/N'^"-), the smallest eigenvalue 
of a A can be as small as 0(1/A^ ") where the constant depends on the smallest eigenvalue of A. For 
a given accuracy e, it is not difficult to show that this algorithm converges after 0(log A^-|-log(l/e)) 
iterations. 

Square root and inverse square root Put C = A^''^ and D = A"^'"^ where A is symmetric 
positive definite. Here, we again choose a constant a such that a|a(x,^)| <C 1 for ^ G {—N,N)'^. 
This also implies that a ~ 0{1/N'^"^). In the following iteration, we first use the Schulz-Higham 
iteration to compute the square root and the inverse square root of aA and then scale them 
appropriately to obtain C and D. 

• Yq = aA and Zq = I. 

• For A; = 0, 1,2, . . ., repeat Ife+i = ^1^^(31 — Z^Yk) and Z^j^i = ^(3/ — ZkYk)Zk until conver- 
gence. 
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• Set C = a'^l'^Yk and D = a^l'^Zk. 

We refer to [23] for a detailed discussion of this iteration. 

In a similar way to the iteration used for computing the inverse, the Schulz-Higham iteration 
is similar to the iteration for computing the inverse in that it uses only additions and products. 
Therefore, all of the computation can be performed via discrete symbol calculus. A similar analysis 
show that, for any fixed accuracy e, the number of iterations required by the Schulz-Higham 
iteration is of order 0(log N + log(l/e)) as well. 

Exponential C = e" . In general, the exponential of an elliptic pseudodifferential operator is not 
necessarily a pseudodifferential operator itself. However, if the data is restricted to ^ G {—N,N)'^ 
and a = 0{1/N '^), the exponential operator behaves almost like a pseudodifferential operator 
in this range of frequencies^ In Section 4.4 we will give an example where such an exponential 



operator plays an important role. 

We construct C using the following "scaling-and-squaring" steps [30] : 



• Pick 6 sufficient small so that a/6 = 2 for an integer K. 

• Construct an approximation Yq for e . One possible choice is the 4th order Taylor expansion: 
Yq = I + SA + ^^ — h ^^1^ — h ^ ^/ . Since S is sufficient small, Yq is quite accurate. 

• For k = 0,1,2, . . . ,K — 1, repeat Ifc+i = Ifcifc- 

• SetC = Yk. 

This iteration for computing the exponential again uses only the addition and product oper- 
ations and, therefore, all the steps can be carried out at the symbol level using discrete symbol 
calculus. The number of steps K is usually quite small, as the constant a itself is of order 0{1/N'^°-). 

Moyal transform Pseudodifferential operators are sometimes defined by means of their Weyl 
symbol aw-, as 

when .^ G Z , otherwise if ^ E M , replace the sum over ^ by an integral. It is a more symmetric 
formulation that may be preferred in some contexts. The other, usual formulation we have used 
throughout this paper is called the Kohn-Nirenberg correspondence. The relationship between the 
two methods of "quantization", i.e., passing from a symbol to an operator, is the so-called Moyal 
transform. The book [I8^ gives the recipe: 

aw[x,i) = {Ma){x,0 = 2" J] j e^^'^--y>^^-^\{y,n)dy, 

and conversely 

a{x,i) = {M~'aw){x,0 = 2" J^ J e-^-'^^-y>^^-^^a{y,v) dy. 



*Note that another case in which the exponential remains pseudodifTerential is when the spectrum of A is real and 
negative, regardless of the size of a. 



19 



These operations are algorithmically very similar to transposition. It is interesting to notice 
that transposition is a mere conjugation in the Weyl domain: a* = M~^(Ma). We also have the 
curious property that 

M^{p, q) = e-^*P«a(p, q) 

where the hat denotes Fourier transform in both variables. 

Applying the operator The last operation that we discuss is how to apply the operator to a 
given input function. Suppose u{x) is sampled on a grid x = {pi/N,p2/N) with < pi,P2 < N. 
Our goal is to compute {Au){x) on the same grid. Using the definition of the pseudo-differential 
symbol and the expansion of a(x,^), we have 

A V 5 

Therefore, a straightforward yet efficient way to compute Au is 

• For each A G {-B^,B^)'^, sample Va(0 for C ^ [-N/2,N/2)^. 

• For each A G {-B^,B^j,f, form the product Va(0(0''"^(0 for ^ ^ [-N/2,N/2f. 

• For each A G { — Bx, B^)"^, apply the fast Fourier transform to the result of the previous step. 

• For each A G {—Bx,Bx)'^, multiply the result of the previous step with ex{x). Finally, their 
sum gives {Au){x). 

Let us estimate the complexity of this procedure. For each fixed A, the number of operations 
is dominated by the complexity of the fast Fourier transform, which is 0(A^log A^). Since there is 
only a constant number of values for A G {—B^, B^)'^, the overall complexity is also 0(A^log A^). 

In many cases, we need to calculate {Au){x) for many different functions u{x). Though the above 
procedure is quite efficient, we can further reduce the number of the Fourier transforms required. 
The idea is to exploit the possible redundancy between the functions /ia,A(C) for different A. We first 
use a rank-reduction procedure, such as QR factorization or singular value decomposition (SVD), 
to obtain a low-rank approximation 

T 

t=i 

where the number of terms T is often much smaller than the possible values of j. We can then 
write 

T 



t=i 




The new version of applying (Au){x) then takes two steps. In the preprocessing step, we compute 
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• For each A G {-B^,B^f, sample ha,x{0 for ^ G [-N/2,N/2f. 

• Construct the factorization /ia,A(0 ^ Z]t=i ^At^t(0- 

• For each t, compute the function '^^ e\[x)u\t. 

In the evaluation step, one carries out the following steps for an input function u{x) 

• For each t, compute vt{C){C)'^""^{0- 

• For each t, perform fast Fourier transform to the result of the previous step. 

• For each t, multiply the result with '^^e\{x)uxt. Their sum gives {Au){x). 

4 Applications and Numerical Results 

In this section, we provide several numerical examples to demonstrate the effectiveness of the 
discrete symbol calculus. In these numerical experiments, we use the hierarchical spline version of 
the discrete symbol calculus. Our implementation is written in Matlab and all the computational 
results are obtained on a desktop computer with 2.8GHz CPU. 

4.1 Basic operations 

We first study the performance of the basic operations described in Section |3j In the following 
tests, we set R = 6, L = 6, and N = R x 3^ = 4374. The number of samples in O is equal to 677. 
We consider the elliptic operator 

Au := (/ - div(a(x)V))u. 

Example 1. The coefficient a{x) is a simple sinusoid function given in Figure El We use the 
discrete symbol calculus to compute the operators C = AA, C = A~^, and C = A^''^. Table [ij 
summarizes the running time, the number of iterations, and the accuracy of these operations. We 
estimate the accuracy using random noise as test functions. For a given test function /, the errors 
are computed using the following quantities: 

• For C = AA, we use "YA(A/t|f" - 

• For C = ^-\ we use ^^^^p^- 

. For C = A^l\ we use ^^^^f^^- 





tj iterations Time Accuracy 


C = AA 
C = A-^ 

C = A^l^ 


3.66e+00 1.92e-05 
17 1.13e+02 2.34e-04 
27 4.96e+02 4.01e-05 



Table 1: The results of Example 1. 
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Figure 4: The coefficient a{x) of Example 1. 
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Figure 5: The coefficient a{x) of Example 2. 





tj iterations 


Time 


Accuracy 


C = AA 


- 


3.66e+00 


1.73e-05 


C = A^^ 


16 


1.05e+02 


6.54e-04 


C = A^l^ 


27 


4.96e+02 


8.26e-05 



Table 2: The results of Example 2. 

Example 2. In this example, we set a{x) to be a random bandlimited function (see Figure |5]). 
We report the running time, the number of iterations, and the accuracy for each operation in Table 

m 

Tables [T] and [2] show that the number of iterations for the inverse and square root operator 
remain almost independent of the function a(x,^). Our algorithms produce good accuracy with a 
small number of sampling points in both x and ^. Although we have A'^ = 4374 in these examples, 
we can increase the value of N easily either by adding several extra levels in the hierarchical 
spline construction or by adding a couple more radial quadrature points in the rational Chebyshev 
polynomial construction. For both of these approaches, the running time and iteration count 
increase only slightly: it is possible to show that they depend on A in a logarithmic way. 



4.2 Preconditioner 

As we mentioned in the Introduction, an important application of the discrete symbol calculus is 
to precondition the inhomogeneous Helmholtz equation: 



Lu :- 



-A 






f 
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where the sound speed c(x) is smooth and periodic in x. We consider the solution of the precon- 
ditioned system 

M-^Lu = M-^f 



with the so-called complex-shifted Laplace preconditioner 

,2 



Ml := -A + 



id 



of which we consider two variants, 

,2 



c^[x) 



and M2 := -A + {l+i) ■ 



UJ 



c^{x) 



For each preconditioner Mj with j = 1,2, we use the discrete symbol calculus to compute 
the symbol of M~ . In our case, applying M~ requires at most four fast Fourier transforms. 
Furthermore, since M-~ only serves as a preconditioner, we do not need to be very accurate about 
applying M~ . This allows us to further reduce the number of terms in the expansion of the symbol 
oiM-\ 

Example 3. The sound speed c{x) of this example is given in Figure pj We perform the test on 
different combination of ui and N with oo/N fixed. For both the unconditioned and conditioned 
systems, we use BICGSTAB and set the relative error to be 10^^. The numerical results are given 
in Table |3j For each test, we report the number of iterations and, in parenthesis, the running time, 
both for the unconditioned system and the preconditioned system with Mi and M2. 




Figure 6: The sound speed c(x) of Example 3. 



(u;/27r, N) 


Unconditioned Mi M2 


(4,64) 


2.24e+03 (8.40e+00) 8.55e+01 (6.40e-01) 5.70e+01 (5.10e-01) 


(8,128) 


5.18e+03 (6.79e+01) 1.50e+02 (4.16e+00) 8.85e+01 (2.46e+00) 


(16,256) 


1.04e+04 (6.50e+02) 4.98e+02 (6.79e+01) 3.54e+02 (4.82e+01) 


(32,512) 


9.00e+02 (6.41e+02) 3.06e+02 (2.20e+02) 



Table 3: The results of Example 3. For each test, we list the number of iterations and the running 
time (in parenthesis). 



Example 4. In this example, the sound speed c(x) (shown in Figure Q is a Gaussian waveguide. 
We perform the similar tests and the numerical results are summarized in Table |4j 

In each of these two examples, we are able to use only 2 to 3 terms in the expansion of the 



symbol of M- . The results in Tables 
the number of iterations by a factor of ! 



and 4 show that the preconditioners Mi and M2 reduce 
to 50, and the running time by a factor of 10 to 25. We 
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Figure 7: The sound speed c(x) of Example 4. 



(w/27r,iV) 


Unconditioned Mi M2 


(4,64) 


3.46e+03 (1.30e+01) 6.75e+01 (5.00e-01) 4.25e+01 (3.20e-01) 


(8,128) 


1.06e+04 (1.39e+02) 2.10e+02 (5.80e+00) 1.16e+02 (3.19e+00) 


(16,256) 


3.51e+04 (1.93e+03) 1.56e+03 (2.16e+02) 6.81e+02 (9.56e+01) 


(32,512) 


1.55e+03 (1.12e+03) 6.46e+02 (4.63e+02) 



Table 4: The results of Example 4. For each test, we list the number of iterations and the running 
time (in parenthesis). 



also observe that the preconditioner M2 outperforms Mi by a factor of 2, in line with observations 
in [15, where the complex constant appearing in front of the up' l(?{x) term in Mi and M^ was 
optimized. 

Let us also note that we only consider the complex-shifted Laplace preconditioner in isolation, 
without implementing any additional deflation technique. Those seem to be very important in 
practice [T5] . 



4.3 Polarization of wave operator 

Another application of the discrete symbol calculus is to "polarize" the initial condition of linear 
hyperbolic systems. We consider the second order wave equation with variable coefficients 

utt — div(a(x)Vn) = 
u(0, x) = uo(x) 
_nt(0,x) = uiix) 

with the extra condition J ui{x)dx = 0. Since the operator L := — div(a;(2;)V) is symmetric positive 
definite, its square root P := L^'^ is well defined. We can use P to factorize the equation into 

(9i + iP)(9i-iP)n = 0. 

The solution n(t, x) can be represented as 

u{t,x) = e**'^n+(x) + e"**^n„(a;) 

where the polarized components u+(x) and U-{x) of the initial condition are given by 



no + [iP) ^ui 



and U- 
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no - {iP) ^ui 



We first use the discrete symbol calculus to compute the operator P ^. Once P ^ is available, the 
computation of u+ and U- is straightforward. 

Example 5. The coefficient a{x) in this example is shown in Figure p] (a). The initial condition 
is set to be a plane wave solution of the unit sound speed: 



uo{x) 



JI-Kikx 



and ui{x) = —2'Ki\k\e 



2TTikx 



where A: is a fixed wave number. If a{x) were equal to 1 everywhere, this initial condition itself 
would be polarized and the component u^{x) would be zero. However, due to the inhomogeneity 
in a{x), we expect both u_|_ and n_ to be non-trivial after the polarization. The real part of u^{x) 
is plotted in Figure [S] (b) . We notice that the amplitude u+(x) scales with the difference between 
the coefficient a{x) and 1. This is compatible with the asymptotic analysis of the operator P for 
large wave number. The figure of u_(x) is omitted as visually it is close to uq{x). 




Figure 8: Example 5. The real part of the polarized component u+ = {uq + {iP)^^ui)/2. Notice 
that the amplitude of u^(x) scales with the quantity a{x) — 1. U- = (uq — {iP)~^ui)/2 is omitted 
since visually it is close to uq. 



Example 6. The coefficient a(x) here is a random bandlimited function shown in Figure^ (a). 
The initial conditions are the same as the ones used in Example 5. The real part of the polarized 
component u+{x) is shown in Figures |9] (b) . Again, we see that the dependence of the amplitude 
of ti+(x) on the difference between a{x) and 1. 

4.4 Seismic depth migration 

The setup is the same as in the Introduction: consider the Helmholtz equation 



+ A_L + 



CO 



c'^{x) 



u 







(13) 



for z > 0. The transverse variables are either x G [0, 1] in the ID case, or x G [0, 1]^ in the 2D case. 
(Our notations support both cases.) Given the wave field u{x, 0) at z = 0, we want to compute the 
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Figure 9: Example 6. The real part of the polarized component «+ = {uq + {iP)~^ui)/2. Notice 
that the amplitude of Ujf-{x) scales with the quantity q(x) — 1. u_ = (uq — {iP)~^ui)/2 is omitted 
since visually it is close to uq. 



wavefield for z > 0. For simplicity, we consider periodic boundary conditions in x or {x,y), and no 



right-hand side in (13). 



As mentioned earlier, we wish to solve the corresponding SSR equation 

|-B(.)),, = 0, 



(14) 



where B{z) is a regularized square root of — A_l — uP' j(?{x^z\ Call ^ the variable(s) dual to x. 
The locus where the symbol 47r^|^p — uP'/c^{x,z) is zero is called the characteristic set of that 
symbol; it poses well-known difficulties for taking the square root. To make the symbol elliptic 
(here, negative) we simply introduce 



a{z-x,i)=g[A^^\i\\ 



1 



UJ 



UJ 



2c^[x,z)) c^(x,z)^ 

where g{x, M) is a smooth version of the function min(x, M). Call b{z; x, ^) the symbol-square-root 
of a{z; X, $,), and B{z) = b{z; x, iVx) the resulting operator. A large-frequency cutoff now needs to 
be taken to correct for the errors introduced in modifying the symbol as above. Consider a function 
x{x) equal to 1 in (— cxo, —2], and that tapers off in a C°° fashion to zero inside [— l,oo). We can 
now consider x{b{z;x,S^)) as the symbol of a smooth "directional" cutoff, defining an operator 
X = xib{z; ,x, —iVx)) in the standard manner. The operator B(z) should then be modified as 

XB{z)X. 

At the level of symbols, this is of course (x(^))tt^tt(x(&)) ^^^ should be realized using the composition 
routine of discrete symbol calculus. 

Once this modified square root has been obtained, it can be used to solve the SSR equation. It 
is easy to check that, formally, the operator mapping n(x,0) to u{x,z) can be written as 

E{z) = exp / B{s) ds. 
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If B{s) were to make sense, this formula would be exact. Instead, we substitute XB(s)X for B{s), 
and compute E{z) using discrete symbol calculus. We intend for z to be small, i.e., comparable 
to the wavelength of the field u(x,0), in order to satisfy a CFL-type condition. With this type 
of restriction on z, the symbol of E(z) remains sufficiently smooth for the DSC algorithm to be 
efficiently the integral over s can be discretized by a quadrature over a few points, and the operator 
exponential can be realized by scaling-and-squaring as explained earlier. 

The effect of the cutoffs X is to smoothly remove 1) turning rays, i.e, waves that would tend to 
travel in the horizontal direction or even overturn, and 2) evanescent waves, i.e., waves that decay 
exponentially in z away from z = 0. This is why X is called a directional cutoff. It is important 
to surround B with two cutoffs to prevent the operator exponential from introducing energy near 
the characteristic set of the generating symbol 47r^|^p — to'^ / c^ (x , z) . This precaution would be 
hard to realize without an accurate way of computing compositions (twisted product). Note that 
the problem of controlling the frequency leaking while taking an operator exponential was already 
addressed by Chris Stolk in |j37j, and that our approach provides another, clean solution. 

We obtain the following numerical examples. 

Example 7. Let us start by considering the ID case. The sound speed c(x) in this example is a 



Gaussian waveguide (see Figure 10 (a)). We set to to be 100 • 2tt in this case. 

We perform two tests in this example. In the first test, we select the boundary condition u{x, 0) 
to be equal to one. This corresponds to the case of a plane wave entering the waveguide. The 



solution of (14) is shown in Figure 10 (b). As z grows, the wave front starts to deform and the 
caustics appears at x = 1/2 when the sound speed c(x) is minimum. 

In the second test of this example, we choose the boundary condition n(x, 0) to be a Gaussian 
wave packet localized at x = 1/2. The wave packet enters the wave guide with an incident angle 



about 45 degrees. The solution is shown in Figure 10 (c). Even though the wave packet deforms its 



shape as it travels down the wave guide, it remains localized. Notice that the packet bounces back 
and forth at the regions with large sound speed c(x), which is the result predicted by geometric 
optics in the high frequency regime. 

Example 8. Let us now consider the 2D case. The sound speed used here is a two dimensional 
Gaussian waveguide (see Figure [IT] (a)). We again perform two different tests. In the first test, 
the boundary condition u{x,y,0) is equal to a constant. The solution at the cross section y = 1/2 



is shown in Figure 11 (b). In the second test, we choose the boundary condition to be a Gaussian 



wave packet with oscillation in the x direction. The packet enters the waveguide with an incident 



angle of 45 degrees. The solution at the cross section y = 1/2 is shown in Figure 11 (c). Both of 
these results are similar to the ones of the one dimensional case. 

5 Discussion 

5.1 Other domains and boundary conditions 

An interesting question is what form discrete symbol calculus should take when other boundary 
conditions than periodic are considered, or on more general domains than a square. 

One can speculate that the discrete Sine transform (DST) should be used as ex for Dirichlet 
boundary conditions on a rectangle, or the discrete Cosine transform (DCT) for Neumann on a 



^For larger E{z) would be a Fourier integral operator, and a phase would be needed in addition to a symbol. We 
leave this to a future project. 
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(b) 



(c) 



Figure 10: Example 7. (a) sound speed c(x). (b) the solution when the boundary condition u{x^ 0) 
is a constant, (c) the solution when the boundary condition m(2;,0) is a wave packet. 

rectangle. Whatever choice is made for e\ should dictate the definition of the corresponding fre- 
quency variable £,. A more robust approach could be to use spectral elements for more complicated 
domains, where the spectral domain would be defined by Chebyshev expansions. One may also 
imagine expansions in prolate spheroidal wavefunctions. Regardless of the type of expansions cho- 
sen, the theory of pseudodifferential operators on bounded domains is a difficult topic that will 
need to be understood. 

Another interesting problem is that of designing absorbing boundary conditions in variable 
media. We hope that the ideas of symbol expansions will provide new insights for this question. 

5.2 Other equations 

Symbol-based methods may help solve other equations than elliptic PDE. The heat equation in 
variable media comes to mind: its fundamental solution has a nice pseudodifferential smoothing 
form that can be computed via scaling-and-squaring. 

A more challenging example are hyperbolic systems in variable, smooth media. The time- 
dependent Green's function of such systems is not a pseudodifferential operator, but rather a 
Fourier integral operator (FIO), where e^'^^^'^ a{x , ^) needs to be replaced by e*(^'^)a(x,^). We 
regard the extension of discrete symbol calculus to handle such phases a very interesting problem, 
see |10| [TT] for preliminary results on fast application of FIO. 
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Figure 11: Example 8. (a) sound speed c(x). (b) the solution at the cross section y = 1/2 when 
the boundary condition u{x,y,0) is a constant, (c) the solution at the cross section y = 1/2 when 
the boundary condition u{x, y, 0) is a wave packet. 



A Appendix 

Proof of Lemma [1} As previously, write 



-2TTix-\ 



' a{x,^)dx 



for the Fourier series coefficients of a{x, ^) in x. Then we can express ^ as 

iAf)ix) = Y^ e2--€ ^ e2--^aA(e)/(0- 



We seek to interchange the two sums. Since a{x, ^) is differentiable d' times, we have 



d\ 



(l + |2vrAr)aA(e) 



^—2nix-X 



[0,1]'* 



{l + {-A,f/^)a{x,Odx, 



hence \ax{0\ < (1 + |2vrA|'^')-i||(l + (-A^)'^'/2)q(2,^^)||^^_ The exponent d' is chosen so that dx{0 
is absolutely summable in A G Z . If in addition we assume / G £i(Z ), then we can apply Fubini's 
theorem and write 

{Af){x) = ^ A' fix), 
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where A f{x) = e^'^*^ (M^^(-^)/)(x), and Mg is the operator of multipUcation by g on the ^ side. 
By Plancherel, we have 

PVIIl^ = \\M^,(of h2 < sup\axm ■ II/IIl2. 
Therefore, by the triangle inequahty, 

\\AfU2 < Y. wA'fu. 



<Y.{1 + \2^\f)-^ ■ sup 1(1 + {-A,f/Mx,0\ ■ II/IIl^ 



As we have seen, the sum over A converges. This proves the theorem when / is sufficiently smooth; 
a classical density argument shows that the same conclusion holds for all / G ^^([0, l]'^). 

D 
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